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abstract. The diffusion/relaxation behavior of polarized spins of pore filling fluid, as often 
probed by NMR relaxometry, is widely used to extract information on the pore- geometry. Such 
information is further interpreted as an indicator of the key transport property of the formation 
in the oil industry. As the importance of reservoirs with complex pore geometry grows, so does 
the need for deeper understanding of how these properties are inter-related. Numerical mod- 
eling of relevant physical processes using a known pore geometry promises to be an effective 
tool in such endeavor. Using a suite of numerical techniques based on random-walk (RW) and 
Lattice-Boltzmann (LB) algorithms, we compare sandstone and carbonate pore geometries in 
their impact on NMR and flow properties. For NMR relaxometry, both laboratory measurement 
$-H and simulation were done on the same source to address some of the long-standing issues in 

its borehole applications. Through a series of "numerical experiments " in which the interfacial 
relaxation properties of the pore matrix is varied systematically, we study the effect of a variable 
surface relaxivity while fully incorporating the complexity of the pore geometry. From combined 
RWandLB simulations, we also obtain diffusion-convection propagator and compare the result 
with experimental and network-simulation counterparts. 
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1. Introduction 

Despite the long history of research on porous media, there remain open issues 
which critically affect various industrial endeavor such as oil/gas exploration, CO2 se- 
questration, water management, storage and migration of toxic waste. A wide range of 
scales permeates through these disciplines, but pore-scale physical processes remain 
their common denominator. In this work, we report our recent effort on the pore-level 
modeling based on micro-tomograms in the context of a borehole application. 




Figure 1 - Pore-grain interface morphology for a typical sandstone (top) and a carbo- 
nate rock(bottom panel). Also shown are the radial pore-to-pore auto-correlation func- 
tion g p (r) at various porosity values (top :three distinct sandstones ; bottom different 
thresholding of the original carbonate tomogram). The insets show scaled correlation 
functions h p (r) which largely collapse into a universal form. 

Several techniques (e.g. | AUZ 96 ]) utilizing detailed 3D pore geometry have rea- 
ched maturity in recent decade thanks to the affordable computing resource, imaging 
techniques, and parallelized simulations. These numerical results are in good standing 
for a class of porous media such as bead packs and sandstones. Here, we focus on 
aspects of more challenging situations involving carbonates in the oil field. 
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Figure 2 - Low-field experimental data (solid zagged curve) and the best simulation 
match (solid curve) with p = 500/im/s. Simplistic model results for NMR response 
at various strengths of surface relaxivity p and pore coarsening (colored broken lines, 
four curves for each value of p) are also shown for comparison. The inset shows the 
porosity and S/V p distribution at the five stages of coarse graining. Shown on the right 
panel are the local magnetization density evolution at various stages of the simulation 
for p = 500/im/s which yields results matching the experiment. 



2. Pore geometry and open issues for the carbonate 



Figure [T] shows two contrasting images of the pore-grain interface for a Fontai- 
nebleu sandstone and a carbonate rock. While quasi-periodicity is clearly visible in 
the former, the latter displays pronounced heterogeneity. This is quantified in the 
right column where the radial pore-to-pore autocorrelation function g(\r2 — ri|) =< 
(/>(ri)</>(r 2 ) > (0(r) = 1 in pore, in grains) is plotted for different values of poro- 
sity. Note that in both cases, the curves collapses to a generic form (insets) h p (r) = 
(g(r) — (j) 2 )/((j)(l — (/))) with <p =< <p(r) >, average porosity. Quasi-periodicity of 
the sandstone is manifest in the form of mild bumps in g(r), reminiscent of the classic 
density correlation in a simple liquid, while the carbonate sample displays a monoto- 
nic, featureless profile, which actually arises from preponderance of many competing 
length scales. In many carbonate rocks, the complexity extends even further than sug- 
gested in the figure, as there exist structures on finer scales beyond the tomographic 
resolution, as well as on scales much beyond the typical sample sizes. Quantitative 
elucidation of both these aspects remain an open challenge which we aim to address 
via further extending the steps described in the following. 
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Figure 3 - The internal magnetic field of a carbonate packstone numerically obtained 
for the entire tomogram (1.5 3 cm 3 ) of the carbonate rock. Shown on the right are the 
slice-cut views of the field and the probability distribution of its strength for a sub- 
block of 512 3 voxels at the center rendered separately for the pore (filled black) and 
grain (gray) space. The curves represent the best Lorentzian (broken lines) | CHE 05 ] 
and Gaussian (solid lines) fits. The latter works better for the pore portion, while nei- 
ther of the methods successfully fits the grain portion. 



3. Numerical NMR relaxometry 

Simulations based on realistic 3D pore allows us to address some of the long- 
standing issues in probing its geometry. One surrounds the validity of the conventional 
mapping between the NMR relaxation spectrum and the pore-size distribution 1KLE 991 
IGRE 071 of carbonate rocks in the possible presence of haphazard heterogeneity in 
their interfacial properties. IRYU 9b 1 The low field NMR data of the carbonate sample 
of Fig [I] (black zagged line in Fig|2]for which bulk relaxation rate of 1/T& = 0.67s -1 
was used and 5 -points moving average was taken), shows the typical multi-exponential 
characteristics. This often invites an interpretation in terms of a broad pore size (more 
precisely, the surface-to-volume (S-V) ratios) distribution. In this scenario, the pore 
space is approximated as a collection of isolated pores of varying sizes, all in the 
so-called fast diffusion limit lBRO 791 . From the tomogram of the same piece, we 
obtain the distribution V of local {Si/Vi} and porosity {Vi} (shown in the inset 
of Figure [2} at various stages of coarse-graining (i labeling sub-blocks of linear di- 
mension L). Attempts to fit experimental data using such recipe neglect the diffu- 
sive coupling between pores as shown by series of curves in the figure with a range 
of controlling surface relaxivity parameter, p, values [ KLE 99 ] (broken lines with 
p = 20, 60, 120 pm/s) and they fail to work : if p is chosen to match the initial slope 
(red curves), it fails to match the experimental data at long times, and vice versa (blue 
lines). This suggests that the model neglecting the extended nature of the pore and the 
heterogeneous diffusive-coupling among its constituents has limited validity in these 
types of rocks. The role of the latter had been previously considered in simple ID 
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models. [RAM 99IIZIE 02l Random-walk based simulation BRYU 0811 on the 3D tomo- 
graphic pore yields an excellent overall agreement (solid black curve) with the choice 
of 500 fim/s for the single parameter. The large p value thus inferred partly accounts 
for the fact that at the resolution of 17 pm per voxel, the digital representation of the 
interface significantly under-estimates the surface area. The effective p, which repre- 
sents the combination of raw p value and the actual surface-area, dictates the long time 
scale dynamics responsible for the good agreement with the entirety of data. 

These numerics can be extended for sophisticated NMR probes[ SON 00 , SON 08 ]. 
The internal field arising from the weak susceptibility contrast between the rock matrix 
and fluid, often a nuisance for NMR probes, offers a way to probe length scales thanks 
to the close geometrical correlation between its spatial profile and the geometry of 
the matrix. I SON 001 Figure [3] shows an example of the internal field profile in the 
carbonate rock calculated using a weak dipole field ansatz, a method verified via direct 
imaging on a pack of cylindrical tubes. |CHO 09 ] We further derive the local field 
gradients that play a critical role in stimulated-echo probes IS ON 001 as well as NMR 
at high fields. II ANA 07NRYU 09a l We find that the field inside the 3D pore space can 
be better approximated by a Gaussian distribution rather than a Lorenzian as reported 
by Chen et a/. [CHE 05 ] A detailed study on the statistics of internal field and their 
gradients inside various types of rocks is under way, as well as their effect on the 
NMR response. 




Figure 4 - Local flow speed through the pore matrix (Fontainebleu sandstone (left) 
and a carbonate packstone (right)), driven from left to right. The color represents the 
local speed and its scheme was chosen to enhance its features. In both cases, the voxels 
with less than 2 % of the maximum speed are omitted from view for enhanced views. 



4. Diffusion-Flow propagators 

One of the objectives of our numerical modeling is to improve the link between the 
static pore geometry (as inferred in situ from borehole measurements) and its trans- 
port, since the latter ultimately controls viability of a hydrocarbon reservoir. The stark 
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contrast between the sandstone and the carbonate pore geometry (Fig[T]) underlies 
the distinct way how the flow is distributed. Figure [4] shows the local flow speed 
distribution for the sub- volume of the rocks obtained from the Lattice Boltzmann 
simulations. 1SUC Oil Clearly, the inhomogeneous flux distribution and extreme tor- 
tuosity as apparent in the carbonate sample hints at why any framework based soley 
on parallel channels of effective hydraulic radii should fail. 



Subvolume of network 




Snapshots from LB+RW simulations 




Figure 5 - The first panel shows the carbonate rock as used in the network simulation. 
The second panel shows the random-walk-LB based propagator simulation on the sub 
volume of the same rock at two different stages of mean displacement. 

The flow propagator BGEN 831 ISTA 841 IARQ 8411 which stands for the probabi- 
lity distribution of tagged particle displacements at two different wait-times (as indi- 
cated by the average displacement Co), P((, Co) provides further details on the in- 
terplay between the geometrical restriction and the flow, and can be measured by 
NMRIHU L9T1 ILEB 961 ISCH 0411 . We simulate its process using two complemen- 
tary techniques : one based on network reduction, the other through a combination 
of random walk and lattice-Boltzmann. (See Figure [5} The former addresses much 
bigger volumes, while the latter incorporates detailed diffusion/fluid dynamics at fi- 
ner length scales. The top-left panel of Figure [6] shows the experimental data on a 
clean dolomite rock. (The rock has a much simpler pore geometry than the carbonate 
rock used in Figs. 1-4.) From the tomogram of the same source rock, we first ran a 
network-model simulation. |ZH A 09 ] The result (lower-left panels) under-represents 
the sharp peak near the origin present in the experimental data. An ad-hoc time de- 
lay factor was introduced, which to a certain degree enhances the weight near zero 
displacement lZHA 091 . Questions arise as to what degree one requires the presence 
of sub-micron-pores which lie beyond the resolution of the tomogram, and how large- 
scale heterogeneity and a finite p affect P(C, Co) - To clarify these, we developed an al- 
ternative method by allowing the random walkers to ride on the background flow from 
an LB run. (the second panel of Fig|5]) Preliminary results (the 2nd and 3rd columns 
for p = and p = 40pm/s) based on the resolution ~ 3pm /vox seem to capture 
the most salient features of the data. This improved agreement provides a valuable 
insight on how one should attribute such experimental features to the heterogeneity at 
extreme length scales. More quantitative inquiry into this issue is in progress. 

The authors would like to thank Drs. J. Dunsmuir (Brookhaven NL), J. Goebbels 
(BAM) and M. Knacks tedt (ANU) for microtomograms of various rocks. 
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Figure 6 - Top panel on the first column shows experimental P(C> Co) at two dif- 
ferent stages. The lower panels show the result of a network model simulation with 
and without the ad hoc time dealy. The middle column shows the combined RW-LB 
simulated propagator with p = at three different stages (parametrized by the ave- 
rage displacement Co)- The third column shows the same but with p = 40/im/s. The 
actual p, from more recent NMR relaxometry, is estimated to be about 20/im/s. Also 
indicated for reference is the diffusion length a t = \JlDt for each stage. 



5. Bibliographic 

[ANA 07] Anand V., Hirasaki G., « Paramagnetic relaxation in sandstones : distingui- 
shing Tl and T2 dependence on surface relaxation, internal gradients and dependence on 
echo spacing », SPWLA 48th Annual Logging Symposium, vol. 446359V, 2007. 

[ARO 84] ARONOVITZ J., Nelson D., « Anomalous diffusion in steady fluid flow through a 
porous medium », Pkys. Rev. A, vol. 30, n° 4, 1984, page 1948. 

[AUZ96] Auzerais R, Dunsmuir J., Ferreol B., Martys N., Olson J., Ramakri- 
SHNAN T., Rothman D., Sch WARTS L., « Transport in sandstone : A study based on 
three dimensional microtomography », Geo. Phys. Lett., vol. 23, 1996, page 705. 

[BRO 79] BROWNSTEIN K., TARR C, « Importance of classical diffusion in NMR studies of 
water in biological cells », Phys. Rev. A, vol. 19, 1979, page 2446. 

[CHE 05] Chen Q., Marble A., Colpitts B., Balcom B., « The internal magnetic field 
distribution, and single exponential magnetic resonance free induction decay, in rocks », 
Journal of Magnetic Resonance, vol. 175, n° 2, 2005, p. 300-308. 

[CHO 09] Cho H., Ryu S., Ackerman J., SONG Y.-Q., « Visualization of inhomogeneous 
local magnetic field gradient due to susceptibility contrast », /. Mag. Res., vol. 198, 2009, 



8 GEOX2010. 



page 88. 

[GEN 83] DE GENNES P., « Hydrodynamic dispersion in unsaturated porous media », 

Journal of Fluid Mechanics Digital Archive, vol. 136, n° 1, 1983, p. 189-200, 
10. 1017/S00221 120830021 16. 

[GRE 07] GREBENKOV D., « NMR survey of reflected Brownian motion », Rev. Mod. Phys., 
vol. 79, 2007, p. 1077-61. 

[HUL 91] HULIN J., Guyon E., Charlaix E., Leroy C, MAGNICO P., « Abnormal diffu- 
sion and dispersion in porous media », Physica Scripta, vol. 1991, n° T35, 1991, page 26. 

[KLE99] KLEINBERG R., « Methods in the physics of porous media », WONG P., Ed., 
Nuclear Magnetic Resonance, vol. 35, Academic Press, 1999. 

[LEB 96] Lebon L., Oger L., Leblond J., Hulin J., Martys N., Schwartz L., « Pul- 
sed gradient NMR measurements and numerical simulation of flow velocity distribution in 
sphere packings », Physics of Fluids Physics of Fluids Phys. Fluids, vol. 8, n° 2, 1996, 
p. 293-301. 

[RAM 99] Ramakrishnan T. S., Schwartz L. M., Fordham E. J., Kenyon W. E., 
Wilkinson D. J., « Forward Models for Nuclear Magnetic Resonance in Carbonate 
Rocks », The Log Analysist, vol. 40, 1999, page 260. 

[RYU 08] Ryu S., « Effects of Spatially Varying Surface Relaxivity and Pore Shape on NMR 
Logging », SPWLA Proceedings of the 49th Annual Logging Symposium, SPWLA, , 2008, 
page 737008 BB. 

[RYU 09a] Ryu S., « Effect of inhomogeneous surface relaxivity, pore geometry and internal 
field gradient on NMR logging : exact and perturbative theories and numerical investiga- 
tions », SPWLA 2009, , 2009, page JJJJ. 

[RYU 09b] Ryu S., Johnson D., « Aspects of diffusive-relaxation dynamics with a non- 
uniform, partially absorbing boundary in general porous media », Phys. Rev. Lett., vol. in 
press, 2009. 

[SCH04] Scheven U., Seland J., Cory D., « NMR propagator measurements on flow 
through a random pack of porous glass beads and how they are affected by dispersion, 
relaxation, and internal field inhomogeneities », Phys. Rev. E, vol. 69, 2004, page 021201. 

[SON 00] Song Y.-Q., Ryu S., Sen P., « Determining multiple length scales in rocks », 
Nature, vol. 406, 2000, page 178. 

[SON 08] Song Y.-Q., Zielinski L., Ryu S., « Two-Dimensional NMR of Diffusion Sys- 
tems », Phys. Rev. Lett., vol. 100, 2008, page 248002. 

[STA 84] Stanley H. E., Coniglio A., « Flow in Porous Media : The 'Backbone' Fractal 
at the Percolation Threshold », Physical Review B, vol. 29, 1984, page 522. 

[SUC 01] SUCCI S., The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford 
University Press, 2001. 

[ZHA 09] Zhao W., Picard G., Leu G., Singer P., « Characterization of single phase flow 
through carbonate rocks : Quantitative comparison of NMR flow propagator measurements 
with a realistic pore network model », Trans. Porous. Media, , 2009. 

[ZIE02] Zielinski L., Song Y.-Q., Ryu S., Sen P., « Characterization of coupled pore 
systems from the diffusion eigenspectrum », /. of Chem. Phys., vol. 117, n° 11, 2002, 
p. 5361-5365. 



